/*
 * TRISpH.cpp
 *
 *  Created on: 16 Aug 2011
 *      Author: Allan
 */

#include "TRISpH.h"

// C++ includes
#include <iomanip>
#include <iostream>
#include <vector>
#include <string>
using namespace std;

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "../../Applications/Equilibrator.h"
#include "../../Assemble.h"
#include "../../AqueousSolution.h"
#include "../../Multiphase.h"
#include "../../Reaction.h"
#include "../../ReactionSystem.h"
#include "PitzerTris.h"

const Multiphase
AssembleTRISMultiphase();

const ReactionSystem
AssembleTRISReactionSystem();

const Equilibrator
AssembleTRISEquilibrator();

void TRISpH(double T, double P)
{
	Equilibrator equilibrator = AssembleTRISEquilibrator();
	
	cout << setw(15) << left << "mTRIS";
	cout << setw(15) << left << "mNaCl";
	cout << setw(15) << left << "pH";
	cout << setw(15) << left << "ln(gTRIS)";
	cout << endl;

	BOOST_FOREACH(double mTRIS, Range(0.0, 0.02, 0.005)) 
	{
		BOOST_FOREACH(double mNaCl, Range(0.0, 5.0, 1.0))
		{
			equilibrator.ClearConditions();
			
			equilibrator.AddCondition("charge-balanced");
			equilibrator.AddCondition("concentration", "H2O", 55.508);
			equilibrator.AddCondition("total-concentration", "Na[+]", mNaCl);
			equilibrator.AddCondition("total-concentration", "Cl[-]", mNaCl);
			equilibrator.AddCondition("total-concentration", "TRIS", mTRIS);
			
			equilibrator.Solve(T, P);
			
			cout << setw(15) << left << mTRIS;
			cout << setw(15) << left << mNaCl;
			cout << setw(15) << left << equilibrator.Get("pH");
			cout << setw(15) << left << std::log(equilibrator.Get("activity-coeff", "TRIS"));
			//cout << setw(15) << left << equilibrator.Get("activity-coeff", "TRIS");
			//cout << setw(15) << left << equilibrator.Get("activity-coeff", "TRISH[+]");
			//cout << setw(15) << left << equilibrator.Get("activity-coeff", "H[+]");
			cout << endl;
		}
		
		cout << endl;
	}
}

const Multiphase AssembleTRISMultiphase()
{
	vector<string> aqueousSpecies = {"H2O", "Na[+]", "Cl[-]", "H[+]", "TRIS", "OH[-]", "NaCl(aq)", "HCl(aq)", "NaOH(aq)", "TRISH[+]"};
	
	AqueousSolution sol(aqueousSpecies);
	
	sol.SetActivityModel("TRIS", PitzerTRIS);
	sol.SetActivityModel("H[+]", PitzerHydron);
	sol.SetActivityModel("TRISH[+]", PitzerTRISH);
	
	Multiphase multiphase({Phase(sol)});
	
	return multiphase;
}

const ReactionSystem AssembleTRISReactionSystem()
{
	Reaction reactionTRISH;
	
	reactionTRISH.SetEquation("TRISH[+]", {"TRIS", "H[+]"}, {1.0, 1.0});
	reactionTRISH.SetEquilibriumConstantData( EquilibriumConstantData(EqulibriumConstantTRISH) );
	
	vector<Reaction> reactions;
	
	reactions.push_back(AssembleReaction("OH[-]"));
	reactions.push_back(AssembleReaction("NaCl(aq)"));
	reactions.push_back(AssembleReaction("HCl(aq)"));
	reactions.push_back(AssembleReaction("NaOH(aq)"));
	reactions.push_back(reactionTRISH);
	
	ReactionSystem eReactions(reactions);
	
	return eReactions;
}

const Equilibrator AssembleTRISEquilibrator()
{
	Multiphase multiphase = AssembleTRISMultiphase();
	
	ReactionSystem eReactions = AssembleTRISReactionSystem();
	
	Equilibrator equilibrator(multiphase, eReactions);
	
	return equilibrator;
}
